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ABSTRACT 

Considering galaxies as self- gravitating systems of many collisionless particles 
allows to use methods of statistical mechanics inferring the distribution function of 
these stellar systems. Actually, the long range nature of the gravitational force con- 
trasts with the underlying assumptions of Boltzmann statistics where the interactions 
among particles are assumed to be short ranged. A particular generalization of the 
classical Boltzmann formalism is available within the nonextensive context of Tsal- 
lis q - statistics, subject to non - additivity of the entropies of sub -systems. Assuming 
stationarity and isotropy in the velocity space, it is possible solving the generalized 
collsionlcss Boltzmann equation to derive the galaxy distribution function and den- 
sity profile. We present a particular set of nonextensive models and investigate their 
dynamical and observable properties. As a test of the viability of this generalized con- 
text, we fit the rotation curve of M33 showing that the proposed approach leads to 
dark matter haloes in excellent agreement with the observed data. 

Key words: galaxies : kinematic and dynamics - dark matter - galaxies : spiral 



1 INTRODUCTION 

Consider a galaxy as an ensemble of many particles evolv- 
ing under the action of their own gravitational potential. 
The enormous number of constituents (~ lC^-lO 11 ) pre- 
vents the solution of the N - body problem both, analyti- 
cally and numerically. Hence, along the history of galaxy 
dynamics studies, different methods have been developed to 
find a way modeling density profiles of galaxies and derive 
their kinematic and observable properties. Notwithstanding 
the fine details, such techniques may be divided into three 
broad classes. 

First, one can rely on an empirical approach starting 
from observations and deducing a theoretical model. For in- 
stance, photometric measurements allow to track the galax- 
ies surface brightness profile, which can be deprojected un- 
der suitable assumptions, yielding the three dimensional lu- 
minosity density. Subsequently, kinematic measurements or 
stellar population synthesis codes allow to determine the ac- 
tual ingredients (e.g., the stellar mass - to - light ratio) lead- 
ing finally to a particular mass model. Typical example s of 
this approach are the exponential model (Freeman 1970) for 



spiral galaxies discs and the PS (|Prugniel fc Simien 19971 ) 



profile for the luminous component of early -type galaxies. 



Based on observations, these methods certainly cannot 
be used to invent models for dark haloes, galaxies are sup- 
posed to be embedded in. This naive consideration has moti- 
vated the development of increasingly high accurate numeri- 
cal N - body simulations, which follow the evolution of many 
particles taking care of the nonlinearities induced by the 
gravitational interactions of each particle with the rest of the 
system. Based on the analysis of the results of these simula- 
tions, a plethora of models for the dark matter component 
has been proposed along the years, typically subject to a 
double power - law profile, i.e. p oc (r /r s )~ a (l + r/r s )~^ 3 ~ a \ 
While there is a general consensus that the density scales as 
r~' ! for large r, there is still an open debate on the value of a, 
adopting a = 1 for the popular NFW model (Navarro et al. 
1996, 1997) or a = -1.5 lMoore et al. I Jl99cj);lGhigna et al. I 
|2000l ) for a steeper and a ~ 0.7 (|Power et al. 1120031 ) for a 
shallower model. It is also possible that the logarithmic den- 
sity slope never attains a finite asymptotic value in the cen- 
tre |Navarro et al. 1 12004| ) th us favouring exponential- like 
mode ls as the Einasto profile (|Einasto II 19651 ; ICardone et al. I 
2003). The still open debate on both, the shape of the den- 
sity profile and its universality along with the disagreement 
regarding the rotation curves of low mass systems (see, e.g., 
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de Blok 2010 and refs. therein) are strong evidences that 
something is still missing or not well understood in the N- 
body simulations, such that the resulting density profiles 
should rather be taken cum grano salts. 

As third approach, we may interpret the observations 
in view of many particle systems, where stars in a galaxy are 
similar to the molecules in a gas. One can therefore rely on 
the methods of statistical mechanics to infer the distribution 
function (DF) of galaxies which, once integrated over the ve- 
locity space, gives the mass density profile. The context of 
standard Boltzmann-Gibbs-Shannon (BGS) statistics along 
with a DF depending on the total energy only leads to the 
conventional isothermal sphere model as equilibrium config- 
uration. However, there is a fundamental difference between 
stars in a galaxy governed by long-range gravitational inter- 
actions and the short range correlations between molecules 
in a gas. Consequently, BGS statistics cannot be applied 
rigorously although it is commonly assumed that the ef- 
fect of this systematic error is negligible. A generalization 
of the BGS statistics has been proposed dealing with sys- 
tems subject to long-range interactions besides other nonlin- 
ear phenomena. Tsallis (1988) first proposed the generalized 
entropy functional of a system as 



p q M +q - 1) 
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(1) 



where k is the Boltzmann constant, pi the probability of 
the i - th microstate and the entropic index q is a real num- 
ber whose deviation from unity parameterizes the departure 
from the standard BGS entropy where 



lim 



i S q = —k^^pi In pi 



(2) 



recovers the classical entropy function. The most distin- 
guishing feature of the generalized nonextensive entropy is 
manifest in the pseudo-additivity as : 

S q (A + B) = S q (A) + S q (B) + (1 - q)k- 1 S q (A)S q (B) (3) 

which again reduces to the standard BGS additivity for 
q = 1. It is worth stressing that the nonlinear term on the 
right hand side accounts for long-range interactions coupling 
two spatially separated systems. Today, nonextensive statis- 
tics in the context of the Tsallis entropy generlizations is 
remarkably successful and widely applied in a variety of dif- 
ferent research fields. 

A first application in galactic dynamics has been per- 
formed by maximizing the Tsallis entropy for fixed to- 
tal mass and energy thus obta ining the analog of poly- 
tropic models for stel l ar sys tems (|Plastino fc Plastino II 1993 ; 
iTaruva fc Sakagami 1 12002 ). Later, Leubner (2005) consid- 
ered the generalized Boltzmann equation and derived den- 
sity profiles for the dark matter haloes in good agre ement 
with N-body simulations |Kronenberger et al. Il200q ). In a 
first approximation Leubner (2005) assumed that the en- 
tropic index q is spatially constant in the system. Indeed, 
since 9 ^ 1 describes the deviations from classical statis- 
tics due to gravitational interactions, it is expected that q 
significantly differs from 1 where the gravitational poten- 
tial is strong, but recovering the classical value q = 1 for 
vanishing potential in the outskirts. Later, a fundamental 
relation combining the spatial dependence of the entropic 
index q with the density and the velocity dispersion was 



provided jDu 1120071 ), without performing a solutions to the 
general problem. Based on this previous analysis, we make 
a step further by investigating the kinematic and observ- 
able properties of a class of models obtained by searching 
for equilibrium configurations in the framework of nonex- 
tensive statistics, assuming spherical symmetry but keeping 
the spatial dependence of the entropic index q along with a 
suitable ansatz for the velocity dispersion required to close 
the system of equations. 

In Sect. 2, we review the derivation of the set of equa- 
tions defining equilibrium configurations in generalized q- 
statistics. After introducing in Sect. 3 a suitable empirical 
form for the velocity dispersion, we also investigate the ba- 
sic kinematics and observable quantities as function of the 
model parameters. As a test case, we fit the rotation curve of 
M33 using the nonextensive approach for the dark halo com- 
ponent. Finally, a summary of results is provided in Sect. 4. 



2 EQUILIBRIUM CONFIGURATIONS 

In the context of theoretical galactic dynamics a system is 
fully characterized by the DF defining the number density 
of stars in the configuration space. Hence, /(r, v, t) is the 
number of particles which, at a given time t, have a position 
in an infinitesimal cube dr centred on r and a velocity in 
a cube dv centred on v. Adopting BGS statistics, the DF 
must be a solution of Boltzmann's equation : 



dt 



+ v 



3/ 
dr 



V*.f£=C(/), 



(4) 



where cj> constitutes the gravitational potential and C(f) is 
the collision term. Assuming a stationary and collisionless 
system, the solution of Eq.(|4]) for spherical symmetry yield s 
the isothermal configuration as ijBinnev fc Tremaine II1987T ) : 



/(r,v) = 



(27TO- 2 ) 3 / 2 



exp 



(5) 



where p 3 denotes the density at the centre, a is the (con- 
stant) velocity dispersion and <j>(r) = —4>{ r ) + <t>o represents 
the relative potential. Integrating over the velocity space 
yields the standard isothermal sphere solution as BGS equi- 
librium configuration. 

For the nonextensive generalization we follow Du (2007) 
defining the q logarithmic and exponential functions as : 
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which reduce to the classical natural logarithm and exponen- 
tial functions in the limit q — > 1 (i.e., when Tsallis stati- 
stics reduces to the BGS case). The generalized Boltzmann 
equation has the same formal content as Q, but the colli- 
sion term is now evaluated in a different way (see Lima et 
al. 2001 and Du 2007 for the full expression). Solving the 
corresponding Lagrangian variational problem for the equi- 
librium configuration, the resulting DF reads : 



f q (e q )^A q [l-(l-q)Pe q ] 1 ^ 1 -^ 
with P being a Lagrange parameter and 

{1 - [1 - (1 - q)/W/2] [1 + (1 ^ q)Pmjj}} 



(6) 



(7) 
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represents the sum of the potential and kinetic energy of 
the system. Note that the total energy is nonextensive ac- 
counting for long-range interactions and correlations in the 
system. Eq.© can be conveniently rewritten in terms of the 
usual variables introducing the velocity dispersion a as : 



/*(r,v) 



B q p a 



(27R7 2 ) 3 / 2 

[1 - (1 - q)v 2 /2o 2 



[l+(l- 9 )V/<T 2 



(8) 



where now a = a(r) is a function of the radial coordinate 
and B q is a normalization constant. For q — > 1, f q reduces 
to the Maxwellian DF and Eq.([8]) denoted the nonexten- 
sive generalization. Integrating over the velocity space the 
density is obtained as : 



P = Ps [l + (1 - q)ip/cr 2 



(9) 



where the potential ip is defined by Poisson's equation for 
spherical symmety : 

1 d 



dr 



>.chp\ 
dr) 



-AirGp 



(10) 



On the other hand, considering Eq.(j8| for an equilibrium 
configuration, Boltzmann's equation Q for a collisionless 
system must be solved leading to a relation between th e 
entropic index q and the velocity dispersion as (|Du 11200^ : 



[do 



GM(r) 



which leads after combining with Eq. (|10[1 to 
aVV + (Vcr) 2 



2ixGp 



(11) 



(12) 



Finally, inserting Eq.© into Eq. (|10[) and eliminating the 
potential ip, yields a second order differential equation for 
the density profile : 



d z p 2 dp q (dp\" 2 fda\ / dp\ A-n 
dr 2 r dr p \dr J a \dr J \dr J G 



4nGp 



q+l 



2„9-l 
Ps 



= 0.(13) 



Eqs. (|12|) and {13} represent a system of two coupled non- 
linear differential equations for the the radial dependence 
of three variables, the density p(r), the velocity dispersion 
a(r) along with the entropic index q(r). When a = const, 
Eqs. (|12|) and (|13[) reduce to those provided by Leubner 
(2005), while, in the limit q — > 1, Eo. (|13p reduces to the 
BGS isothermal sphere condition with a = const such that 
Eq. (|12[) becomes an identity. 



3 SPHERICAL MODELS 

Derived only under the assumption of spherical symmetry, 
Eqs. H12p and (|13p describe the equilibrium configuration of 
a spherically symmetric gravitational system of collisionless 
particles in the framework of Tsallis q - statistics. The system 
is not closed since we have only two equations to determine 
the three quantities p(r),a(r) and q(r). In order to find a 
solution, we have to make an ansatz for one of these quanti- 
ties and check a posteriori whether the corresponding den- 
sity profile is physically reasonable. For instance, we may 
postulate a radial dependence of the entropic index q(r), 



rather unsuitable since this parameter is related to complex 
nonlinear interaction phenomena. Alternatively, it is easier 
to parameterize in a model independent way the form of the 
velocity dispersion looking at the measured profiles in liter- 
ature. We assume the following ex pression for the velocity 
dispersion (|Napolitano et al~ll2010l ) : 



a(r) 



(To 



1 



r + r Q 



(14) 



where (ro, <7o) are scaling quantities and n controls the tran- 
sition from the varying behaviour to the constant asymp- 
totic value. A caveat is in order here. Eo. (|14|) was proposed 
to describe the radial component oy of the velocity disper- 
sion since, when projected along the line of sight (los), it 
assumes a constant <j; os profile in agreement with planetary 
nebulae measurements in NGC 4374. We have here taken 
this same phenomenological expression as an input for the 
cr(r) function entering the DF, interpreting er(r) as a velocity 
dispersion. Actually, such an identification is only a formal 
one, but what physically cr(r) represents is not fully under- 
stood. To understand this point, let us consider the limit 
q = 1 so that the DF reduces to the Maxwellian one. This 
gives rise to the well known isothermal sphere model having 
the remarkable property that the radial and the tangential 
velocity dispersion are equal. In such a case, a(r) turns out 
to be both the total velocity dispersion and the radial one 
(modulo a scaling factor). On the contrary, when q 1, it 
is not clear whether a(r) should be interpreted as the radial 
oy(r) or the total velocity dispersion or as a fully differ- 
ent quantity having the same physical dimensions as oy(r). 
Should one interpret cr(r) in Eq.© as the radial velocity 
dispersion o- r (r), then one could rely on the spherical Jeans 
equation relating ay(r) and p(r) and assume a parametrized 
expression for the anisotropy profile /3(r) = 1 — ae(r)/a r (r). 
Adding this relation to Eqs. (|12|) and (|13p , a closed system 
for the three unknowns p(r), a(r), q(r) is obtained and can 
be solved numerically. Although such an approach is worth 
to be explored, it is nevertheless flawed by the theoretical 
uncertainty on validity of the identification a(r) = oy(r) for 
non Maxwellian DFs and still need an assumption for a func- 
tion, namely the anisotropy profile. We therefore prefer to 
adopt a phenomenological parameterization of a(r) to avoid 
giving a conclusive interpretation of its physical meaning. 
Defining the following dimensionless quantities : 

x = r/r , p = p/p s , cr = a/a Q , 

we can conveniently write the system (|12[) - (|13[) as : 

g +[W .> + *<*,)i2-i(2)' + *£-o.M> 

2/K 



1-9= ■ Sq(x,7]) , 

p 

with : 



/C = 



±nGp s rl 



(16) 



(17) 



1 A similar situation also takes place when one tries to model 
the coarse grained distribution function, defined as p(r)/a 3 (r). 
Depending on whether cr(r) is meant as the total or the radial 
velocity dispersion, one gets different results for the density profile 
and the other quantities of interest. 
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Figure 1. Scaled density profile p = p/p s as a function of the dimensionless radius x = r/ro for models with ao = —0.5 and different 
values of (r), K,). In the left panel, we set K = 15.0 and r/ = 0.1 (short dashed), 1.0 (solid), 2.0 (long dashed), 3.0 (dot dashed). Similarly, 
in the right panel, r\ = 1.0 and we consider models with K = 7.5 (short dashed), 15.0 (solid), 30.0 (long dashed), 45.0 (dot dashed). 



l+[x/{l + x)\n 

S q (x,r,) = ~[aV 2 a + (Va) 2 ] (19) 

{2[x/(l + a)]" - l}y 2 - {1 + [x/(l + x)]"}t, 
x 2 ^(l + x) 2 +*>{1 + [x/(l + a;)]"} 4 

Eq. (|16[) is now a simple algebraic equation for q{rf) to 
be solved trivially and inserted into Eq. (|15[l . This proce- 
dure generates a complex second order nonlinear differential 
equation for the scaled density p, which can be straightfor- 
wardly integrated numerically provided that initial condi- 
tions are given as follows : 

p{Xmin) — 1 ; P \Xmin) — ao/a^mm , 

where x m in is a very small but finite number, the prime 
denotes the derivative with respect to x and ao = 
d\n p/dlnx(x = x m in) is the logarithmic density slope in 
Xmin- Note that one should in principle set x m i n — since, 
by definition, p(0) = 1. However, because of numerical prob- 
lems, we will never solve Eq. (fr5J) up to x min = 0, but 
rather set x m i n — 10~ 5 and assume that the difference 
p{x m in) — p(x = 0) is negligible. As a further remark, we 
note that ao must be zero to assure that the density does 
not diverge at the centre and p a in the DF Q takes a finite 
value everywhere. However, this does not prevent us from 
taking a negative value for ao if it is evaluated in x m i n > 0, 
still imposing ot(x = 0) = 0. While this introduces a non- 
realistic abrupt change in the logarithmic density slope pro- 
file, we nevertheless prefer leaving ao as a free parameter in 
order to explore a larger class of models. Thus, the result- 
ing models are not suitable for x ^ Xm in , but x m i n is small 
enough to avoid any impact on astrophysical application of 
interest. With this caveat in mind, we explore the physical 
properties of the model obtained by solving Eqs. (|15l) - (|16l) as 
a function of the three parameters (r/, K,, ao), while stressing 
that (ro, Co) are only scaling parameters. 

Before discussing the properties of the model thus ob- 
tained, it is worth pointing at an underlying assumption. 
As Eq.Q shows, the DF is isotropic in the velocity space so 



that the velocity distribution function <j>(y) (i.e., the integral 
of the DF over the position coordinates) is the same along 
both the radial and t angential direction. Actually, recent 
numerical simulations |Hansen et al. 1 120061 ; iHansenl |2009| ; 
iKuhlen et al~ll2010l ; [Ling et al. II2OI0F have pointed out at a 
marked difference. To overcome this issue, one should modify 
the dependence of the DF on the velocity which now enters 
only through its magnitude (entering the total energy), but 
not its separated components. In such a case, however, the 
DF would not be the solution of the variational problems 
for the non extensive statistics which should then be refor- 
mulated from the beginning. In this first step analysis, we 
therefore to retain the assumption (f>(v r ) = <l>( v t) warning 
the reader that the resulting model should be used only for 
those subset of dark matter haloes where this ansatz ap- 
proximately holds true. 

3.1 Density profile 

In order to solve numerically Eo. (|15[) . we need to pro- 
vide the parameters (r/, KL, ao) along with the initial con- 
ditions and choose a range for the scaled radius x. We fix 
(x m in, x max ) = (10 -5 , 10) such that the numerical solution 
can be reliably performed avoiding both, the very inner parts 
affected by convergence problems and the very outer regions 
where a significant loss of precision makes the numerical 
solution unstable. As we will see later, the spatial range 
adopted is much wider than the one usually probed by ob- 
servations, wherefore no region of interest is excluded. 

Fig-E shows the scaled density p — p/p 3 as function 
of the dimensionless radius x = r/ro for different choices 
of the model parameter^ (rj, /C, ao). As a general result, we 
find that the density profile decreases everywhere so that the 
model is physically meaningful, but the rate of this scaling 
with x depends on the adopted parameters. 

For fixed values of (rj,IC), the density profile turns out to 
be almost independent of ao, affecting the shape of p(x) only 
in the very inner regions (x « 0.01). This is an expected 
result noting that ao only enters into the initial conditions, 
thus limiting a global influence. As soon as x increases, the 

2 The motivations for the values of (r/, /C, ao) adopted in this and 
the following figures will be discussed later . 
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Figure 2. Same as Fig. [T] but for the 

dominant terms in Bq. (|15[) are driven by (r],IC), making the 
exact value of Qo unimportant. Therefore let us assume ao = 
—0.5 as a compromise between cored and cusped models and 
consider the variation of different quantities of interest as a 
function of (77, IC) only. 

The density profile indeed strongly depends on these 
two parameters. In particular, we note that, except in the 
x « 1 limit, the larger 77 is, the larger is also p(x) for a 
fixed IC, although the different profiles converge to each other 
asymptotically. This effect can be qualitatively explained by 
looking at the assumed velocity dispersion profile in Eg. (1141) . 
While <r(r) always approaches a constant value, its initial de- 
crease is more pronounced for larger 77 values. Since a larger 
(j{x) calls for a larger density, it is therefore clear why p(x) 
has to be an increasing function of 77 in order to reproduce 
the assumed a(f) profile. The impact of 77 is, however, im- 
portant only in a limited radial range and, moreover, there 
seems to be a sort of saturation such that values of 77 << 1 
give rise to models that can be hardly discriminated. 

IC appears as important parameter of the theory and 
is indirectly proportional to the density for fixed x and 77 
values. Understanding the motivation for this behaviour is 
not straightforward. On one hand, Eq. (|17|) shows that the 
central density p„ is larger for larger K, so that one naively 
expect p oc IC should Eq. (|17[) be extrapolated to any radius. 
One can, however, speculate a reason for the inverse scaling 
of p with IC looking at Eq. (|16l) for the entropic index q(x). 
To this end, let us first note that the adopted parametriza- 
tion for cr(r) gives a velocity dispersion which quickly ap- 
proaches a constant value as is the case for the isothermal 
sphere. As well known, this latter model is the equilibrium 
configuration for the BGS statistics characterized by q = 1. 
It is therefore expected that 1 — q is roughly constant and 
small everywhere but in the very inner regions where the 
long range gravitational interactions are maximized. As a 
consequence, the ratio S p / (ICp) must be constant and, since 
S p (x) is roughly constant for x > 0.1 — 0.5, so must K,p(x) 
be if we want that the right hand side of Eo. (|16|) has the 
expected behaviour. As a consequence, the larger IC is, the 
smaller must be p(x). We may therefore argue that the scal- 
ing of p with IC is just a manifestation of Tsallis statistics 
being a moderate correction to the BGS statistics in colli- 
sionless self- gravitating systems. Note that this is related to 
the particular expression adopted for a[r), but we can argue 
that a similar behaviour generally holds true provided the 
DF is not too different from a Gaussian. 



i3 




2 4 6 8 10 



ithmic density slope a = din p/ dlnx. 

While the density profile is always decreasing, the local 
slope of the p vs x relation is not constant. This can be 
clearly appreciated looking at the logarithmic density slope : 

_ d In p d\np 
d\nr dlnx ' 
which is plotted in Fig.[2]for the same choice of model param- 
eters as in Fig. [I] Indeed, a is a strongly varying function of 
x starting from ao, decreasing toa~ —3 for x = x a ~ 1 and 
then slowly approaching the isothermal a = — 2 value (out- 
side the range plotted). It is interesting to note that it now 
77 playing a more important role than IC. In particular, the 
larger is 77, the larger (in absolute value) is the minimum of 
a(x) and the later it is achieved, i.e., the larg More- 
over, for x ^ x a , \a(x)\ is a quickly decreasing function of 
x, while the opposite trend is shown in the x x a range. 
The dependence on IC is reversed with respect to those on 
77 thus mirroring what happens with p(x). We can qualita- 
tively explain the minor role of IC noting that, according 
to our previous interpretation, the model evolves in a way 
keeping ICp nearly constant, whatever the scaling of p with 
x is. Contrary, since 77 determines the transition from the 
decreasing to the flat profile of cr(x), it has to play a similar 
role for p(x), thus explaining why a strongly depends on it. 

The a(x) profile suggests that the density law could 
be approximated by a dou ble power - law mo del such as 
the generalized NFW one (|jing fc Suto 1 12000| ) with p oc 
(r/r 3 )~ a [l + (r/r s )]- (3_a) , or the (a, ,3,7) models (Zhao 
1996, 1997) having p oc (r/r s )"° : [l + (r/r s ) 7 ]- (/3 - a)/7 . Actu- 
ally, all these models present a central divergence, while our 
model is constructed to have p = 1 in x = 10~ J . Moreover, 
there is a single scale radius so that the transition is quite 
smooth, while our a(x) profiles show a quick transition. By 
trial and error, we have found that a better approximation 
is provided by the following profile : 







-1 






1 + 


(sr: 




1 + 


(in 













(20) 



yielding a small core for x « x c , decreasing as x~ ain for 
x c « x « xt and asymptotically scaling as x~ a ° ut in 
the outer regions. The two radii (x c ,x t ) separate the differ- 
ent scaling regions, while a t determines the smoothness of 
the transition. In order to check whether Eq. (|20fl indeed ap- 
proximates well the numerical solution, we have fitted it to a 
large sample of models randomly generating the parameter 
space (ri,IC,ao). A satisfactory fit (with an rms percentage 
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0.001 



10 



Figure 3. Scaled mass M(x)/(ATTp B r^) as a function of x for 
models with (ao,f?) = (—0.5,1.0) and K. = 7.5 (short dashed), 
15.0 (solid), 30.0 (long dashed), 45.0 (dot dashed). 



deviation of the order of 7%) is found for ~ 50% of the sam- 
ple, but one should investigate in more detail which are the 
sets (rj, K., «o) giving rise to realistic models. We therefore ar- 
gue that Eq. (|20p provides a reasonably well fitting formula, 
which we checked in addition by fitting models tailored on 
observed galaxies. For these successful fitted models, median 
and median deviation of the slope parameters are found to 
be: 

a in = 0.22 ± 0.13 , a t = 3.78 ± 1.51 , aw = 3.14 ± 0.10 . 

We also note that about half of the successful fits are subject 
to very large xt values, better approximated by the simpli- 
fied empirical form : 



1+1^ 

X c 



p(x) = 

with a single transition scale and 



(21) 



a in = 5.94 ± 0.28 , a out = 2.41 ± 0.08 . 

Investigating which model is more appropriate and how the 
fitting parameters change with (rj, K,, ao) is outside our aims 
here. We nevertheless note that, whichever model is pre- 
ferred, a(r/ro » 1) < 10/3 in agreement with the upper 
limit found by Hansen et al. (2005) for the case of constant 

q- 



3.2 Mass profile and rotation curve 

Within the restriction of spherical symmetry, the mass pro- 
file of the nonextensive model is trivially evaluated as : 

M(x) = 4tt / r' 2 p(r')dr' = 47rp s rjj / x' 2 p(x')dx . (22) 
Jo Jo 

Fig.[3] shows the scaled mass M(x) — M(x)/(A-Kp s r'l) as a 
function of x for fiducial values of (00,77) = (—0.5, 1.0) and 
different K.) (the same results hold for other values). The 
dependence of M(x) on the model parameters corresponds 
the one of the scaled density, as expected in view of their 
interrelation. Note that M(x) depends much stronger on 77 
than p as could be anticipated, since this quantity is directly 
related to cr(r), while the relation between a and p involves 
an integration which smoothes out the scaling with rj. 

It is worth wondering whether the total mass of the 
model is finite or not. On the one hand, Fig. [3] shows that 



M(x) keeps increasing up to the last point of the radial 
range probed thus suggesting an infinite total mass should 
the trend be extrapolated. This is consistent with a(x) ap- 
proaching the isothermal a — —2 value. Considering the 
shape of the density and logarithmic slope profiles, we ar- 
gue that our model is similar to a cored isothermal sphere in 
the very inner and very outer regions claiming that the total 
mass is indeed infinite. However, Eq. (|15[) cannot be solved 
to radii larger than x max = 10 thus preventing a safe con- 
clusion regarding the above (although reasonable) inference. 
Should this be the case, the model should be restricted to 
a maximum radius, i.e. the viral radius (defined in the way 
that the mean density within R v i r is A v i r pM with A v i r de- 
pending on the cosmological model and pm the mean cosmic 
matter density). As we will see later, x V i r — R v i r /ro is found 
typically well within the upper limit x ma x = 10 implying 
securely that our solution describes the relevant quantities 
over the range of interest for applications. It is worth not- 
ing, however, that such a high radius cutoff is likely to be 
not necessary at all. Indeed, it is still possible that the in- 
finite mass problem (if present at all and not an artifact of 
the numerical solver) is just a consequence of the adopted 
parametrization of o(r). Indeed, Eq. (|14[l has been proposed 
to get an asymptotically flat dispersion at large radii just as 
in the case of the popular isothermal sphere. Indeed, a flat 
cr(r) can alternatively be obtained postulating M(r) a r at 
large radii thus giving rise to a formally infinite mass. It is 
therefore possible that a different choice for cr(r), such as a 
double power -law, leads to a model which is similar to the 
one we are discussing over a large radial range, but presents 
a finite total mass. In order to escape the arbitrariness in 
the choice of a different a(r) profile and to avoid introduc- 
ing further parameters, we will not consider this possibility, 
but warn the reader against concluding that non extensive 
statistics lead to infinite mass models. 

The mass profile serves as input quantity for evaluat- 
ing the rotation curve defined by v 2 (r) — GM(r)/r, i.e. 
v 2 (x) oc M(x)/x. Since M(x) approximately grows linearly 
for x > > 1 , an asymptotically flat rotation curve is achieved 
in agreement with observations in spiral galaxies. However, 
the question arises whether the model we are discussing 
is used to describe the full mass content or the dark halo 
only. For a typical spiral galaxy, the stellar component is 
excellently described by an exponential profile, while the 
nonextensive approach is roughly represented by a double 
power -law, implying that this model is applicable for the 
halo component only. In this case, the entire rotation curve 
is obtained by summing up the stellar and dark matter con- 
tributions where the former likely dominates in the inner re- 
gions. Hence, it is important to tailor the ro value such that 
the regime x >> 1 corresponds in linear units to r >> 
with Rd being the disc half mass radius. In this case a flat 
rotation curve can be generated in the outer regions of the 
disk as it is indeed observed. We will come back to this 
point later when considering the application to an actual 
case, showing how this can be easily achieved. 



3.3 Surface density and projected mass 

Although argued above that our model is not suitable to de- 
scribe the disc of spiral galaxies, it is nevertheless interesting 
to compute the surface mass density, since this quantity en- 
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Figure 4. Scaled surface density E(x)/(p s ro) as a function of 
x for models with (ao,/C) = (—0.5,15.0) and r\ = 0.1 (short 
dashed), 1.0 (solid), 2.0 (long dashed), 3.0 (dot dashed). 



ters in the determination of the lensing potential and hence 
deter mines the magnification and shear induced by the lens 
halo (|Schneider et al. Ill992l ). To this end, we just have to 
replace x = \/ £, 2 + C 2 in the density profile and integrate 
over £, upon assuming that the line of sight coincide with 
the z - axis with (£, () — (R/ro, z/ro). The results are shown 
in Fig.[4]for different choices of r) and (ao,IC) = (—0.5, 15.0) 
as fiducial case. As general feature we note that ~S(x) first 
decreases approximately, changing to a shallower profile as 
x increases. This behavior differs signifficantly from an ex- 
ponential or Sersic (1968) profile wherefore it is not possi- 
ble applying this model for surface brightness fits of spiral 
or early -type galaxies. Moreover, it enforces the previous 
suggestion that this model does not apply for the stellar 
components, but can conveniently be adopted for modelling 
dark haloes. Not surprizingly, concerning the dependence on 
the model parameters we find that the dimensionless surface 
density E(a;) = ~E(x)/(p 3 ro) scales with (rj,IC) in the same 
way as the 3D density p(x). Indeed, this is expected since 
the projection does not invalidate any of the argument pre- 
viously discussed regarding the scaling properties of p. 

A quantity of great interest in lensing applications is 
the deflection angle, which, for a spherically symmetric lens, 
reads a(£) = M pro j(£) /(7r££ crit ), where 



Ah 



(£) = 2tt / Y,(R')R'dR' 



E(0?'d£'(23) 



is the projected mass within the radius R = £ro and the crit- 
ical surface density is defined as £ cr it = c 2 D s /(4ivGDdDds) 
with (Dd,D s ,Dds) being the observer - lens, observer - 
source and lens -source angular diameter distances, respec- 
tively. In case the system observer - lens - source is (nearly) 
perfectly aligned, the image of the lensed galaxy is a ring 
with radius Re = &(Re), referred to as the Einstein ra- 
dius. A measurement of Re then allows to obtain the pro- 
jected mass within Re as M pro j(RE) = irT, cr itR'E provided 
that lens and source redshifts are known and a cosmological 
model has been set to compute angular diameter distances. 

While the dependence of M pro i(C) on the model pa- 
rameters is the same as for the entire mass M(x), it is 
worth stressing that there is a significant difference among 
the various (rj,IC) combinations. Hence, if M pro ,,-(£) is es- 
timated from the Einstein radius in lensed systems, nar- 



row constraints on the model parameters follow. To this 
end, the projected mass for £ > 1 should be measured, 
which actually could be rath er difficult. Consid ering, for 
instance, the SLACS sample (1 Auger et al. I [20091 ') . we find 
(Re I Ref f) ~ 0.6 requiring that ro/i? e // < 1 in order to 
have M pro j(^E) probing the range £ > 1 where the dis- 
crimination is possible. Actually, unless the halo is quite 
concentrated, it is unlikely that ro < Reff leaving us in 
the opposite (£ < 1) regime. However, we might rely on 
weak lensing to constrain M pro j(£) at larger radii or com- 
bine lensing and central velocity disper sion measurements 
(|Cardone et al. ||2009| ; lAuger et al.ll2Q10l ) for discriminating 
better among different (rj, K) values. 



3.4 Aperture velocity dispersion 

In spiral galaxies the rotation curve is easily measured due 
to the presence of large amounts of gas probing the gravita- 
tional potential field. Contrary, early -type galaxies as older 
systems are subject to low gas content where a determina- 
tion of the circular velocity is quite difficult. Moreover, these 
systems are typically dominated by random rather than or- 
dered motions leading to a small v c (r), complicating quan- 
titative measurements even more. On the other hand, the 
velocity dispersion profile may be traced re lying on stars 
(|Cappellari et al. 1 120061 ; iThomas et "all 120091) for the inner - 
regions and, e.g., pla neteary nebulae ( Douglas et al. I [20071 : 
ICoccato et"aT~ll2009l ) in the outskirts. 

In order to determine this quantity, Jeans equations for 
a spherically symmetric systems must be solved first yielding 
the radial velocitjjf). Based on an ansatz for the anisotropy 
pro file the result can be projecte d along the line of sight (see, 
e.g.. lBinnev fc Tremainel (|l987t ) for a detailed description). 
It is worth discussing whether such a procedure still holds 
in our case. Indeed, the Jeans equations follow from the mo- 
ments of the collisionless Boltzmann equation indicating the 
implicit presence of BGS statistics. Actually, Boltzmann's 
equation still holds for nonextensive environments allowing 
to rely on the standard formalism. Hen ce, the los veloc- 
ity d ispersion can still be computed as (|Mamon fc Lokasi 
l2005h : 



GM, 



ffP 



T*I(i?) 



K(x e /£e)p*(€e)Mtot(3: e 



-cte e (24) 



intensity profile, 



where I(R) is the 

(r / R e f f , Rj R e f /), M e ff and pl JJ are the total mass 
and the stellar density, respectively, at the effective radius 
Reff and the hatted quantities are normalized with respect 
to their values at Reff- Finally, Mtot(xe) is the total 
mass, while K(x e /£,e) is a kernel function depending on 
the choice of the anisotropy profile. As hinted at before, 
in our approach, we have implicitly assumed that the 
radial and tangential velocity distribution functions are 
equal. As a consequence, our model is isotropic in the 
velocity space by construction so that /3 = is the only 
possible choice. We take the corresponding expression for 
K(xe/£e) from Appendix B of Mamon & Lokas (2005) and 

3 It is worth stressing that, should one assume that a(r) is yet 
the radial velocity dispersion, this step is redundant and one can 
directly project cr(r) along the line of sight. 
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use a de V aucouleurs (1948) profile for the intensity law 
and a PS (|Prugniel fc Simien I Il997l ) model (thus setting 
p* oc Xe Pn exp[— b n %e ] and n = 4) for the stellar density. 

Since we are not interested in fitting any actual dataset, 
the stellar component is neglected in order to demonstrate 
better the behavior of oi os for our modeQ. Fig.[S]shows some 
profiles for the usual combinations of model parameters and 
two values for the ratio ro/R e ff which adds now to the list 
of quantities cr ios (£ e ) depends on. Since <r; os is related to 
the mass profile, it is not surprising that, for given £ e , its 
dependence on the (77, /C) parameters is the same as that 
of M(x), see the discussion above. We only remark that 
the radial range where the oi os profiles significantly differ 
depends on the adopted value of the ratio ro/Reff- The 
smaller is this ratio, the larger is the radial range where 
discriminating among different models is possible. Moreover, 
o"ios(£) is larger for smaller ro/Reff value as a consequence 
of the interplay between the mass profile and the intensity 
I(R)- Indeed, because of the exponential term in the stellar 
density, the main contributions to the integral in Eq. (|24[) 
comes from the regions within few times R e ff indicating 
that the smaller is ro/Reff, the more the mass term M(x e ) 
contributes and the larger is the los velocity dispersion. 

Measuring ai os (R) is actually possible only if the galaxy 
is quite close providing spectra at different distances from 
the centre. For distant galaxies (for instance at intermedi- 
ate redshift, typical of SLACS lenses), it is only possible 
to measure luminosity weighted velocity dispersion within 
a circular aperture of radius R ap set by the instrumental 
setup. Fig.[S] presents this quantity evaluated as: 

2nJ«^ a los {R)I{R)RdR 

(T a p(Rap/Reff) = " Tr — (25) 

2tt J""* I(R)RdR 

and assuming R ap /R e ff = 1/2, as for the median values in 
SLACS lenses. While the scaling with (77, JC) is the usual one, 
we note that the weighting procedure strongly smoothed out 
the trends, such that discriminating among different cases 
using cr a p only becomes very difficult unless ro ~ R e ff> 
which is quite unlikely. 




Figure 6. Dimensionlcss aperture velocity dispersion as a func- 
tion of the ro / R e f f ratio for a circular aperture of radius R ap = 
R e ff/2. The model parameters are set as in Fig,l5l 



4 TESTING THE MODEL 

The model we have investigated up to now is theoretically 
well founded and motivated by applying Tsallis statistics 
to the case of a collisionless self - gravitating system. How- 
ever, we have to proof that this approach is also reliable 
when applied to a realistic structure, i.e. that this model 
is indeed able to reproduce observations. As a first appli- 
cation, we consider the case of M33, a spiral galaxy at a 
distan ce of 0.7 Mpc with an ex cellently measured rotation 
curve (|Corbelli fc Salucci I I2OO0I ) extending up to 13 times 
the disc half mass radius. 

Before examining the fit to these data, we have to dis- 
cuss preliminary how the galaxy is modelled. As usual, we 
split the total mass in two components, the stellar one and 
the dark matter halo. The surface brightness profile is well 
fitted by an exponential profile so that we assume an in- 
finitesimally thin disc for th e stellar compo nent and com- 
pute the circular velocity as (|Freeman Ill970l ) : 

v 2 d (R) = {2GM d /R d )y 2 [I (y)K (y) - h&K^y)] (26) 



4 Needless to say, this introduces a severe systematic error in the 
inner regions which are typically stars dominated so that the plots 
in Fig.[5]should not be compared to measured profiles. 



where y = R/2R d , In{y) and K n (y) are the modified Bessel 
functions of order n of the first and second kind, respectively, 
and the total disc mass is M d = Y+L d with L d the total 
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Figure 7. Best fit rotation curve superimposed to M33 data. 

luminosity and T* the constant stellar mass - to - light raticQ 
Photometric measurements allow to set Ld = 4.2 x 10 9 L0 
and Rd = 1.2 kpc leaving T* as only parameter. 

Next, adopting the nonextensive context to model the 
dark halo density distribution the question arises whether 
this approach is self consistent. Indeed, the derivation was 
performed with respect to a system of collisionless self- 
gravitating particles, implicitly assuming a single compo- 
nent environment. Actually, we consider now a two compo- 
nent structure and, moreover, the stellar system is described 
by a model which is not the equilibrium configurations of 
a given statistic^]. However, we may assume that our ap- 
proach applies to the proto - galactic systems providing the 
equilibrium configuration for the dark matter particles be- 
fore the onset of baryonic collapse. Then it is still possible 
using the resulting model for the present day halo, provided 
that the disc onset has not significantly modified the start- 
ing configuration of DM particles. This should be a good 
approximation everywhere unless in the very inner regions 
where stars give the main contribution. However, the error 
induced on the total circular velocity should be negligible in 
case that these inner regions are stars dominated. 

With this caveat in mind, we can fit the rotation curve 
of M33 in order to check the viability of the nonextensive 
halo model and constrain the stellar M/L ratio and the halo 
parameters. To this end, we adopt a Bayesian approach and 
maximize the likelihood function £(p) oc exp [— x 2 (p)/2] 
with p = (T*, 77, K., ao, ro, 00) and 



x 2 (p) 



J^abs 

E 

i=l 



'(Ri) 



\Ri,p) 



where v° bs (Ri ) and v^ h (Ri , p) are the observed (with error 
Ei) and theoretically predicted circular velocity at Ri and the 
sum is over the Nobs points. In order to efficiently sample 
the six dimensional parameter space, we use a Monte Carlo 



5 We stress that Y* only refers to the stellar component and 
must not be confused with the total M/L = M(R)/L(R) with 
M(R) = Mdisk(R) + M DM (R) which is not constant. 

6 We note that, actually, application of our formalism to the disc 
is not possible since baryons are subjected not only to gravita- 
tional interactions, but also to complicated phenomena describing 
the formation and evolution of stars. Taking care of all these in- 
gredients in a statistical description is a daunting task that, as 
far as we know, has still to be addressed. 



X 


x B f 


(x) 




68% CL 


95% CL 


T. 


1.21 


1.22 


1.21 


(1.17, 1.26) 


(1.12, 1.31) 


'/ 


5.25 


1.99 


1.20 


(0.33, 2.86) 


(0.13, 10.86) 


K. 


26.27 


20.86 


13.00 


(6.90, 39.93) 


(3.25, 74.22) 


ro 


13.08 


15.06 


13.24 


(8.92, 21.92) 


(6.40, 32.13) 


TO 


76.3 


101.5 


96.1 


(79.9, 128.3) 


(53.7, 146.6) 



Table 1. Constraints on the model parameters from the fit to 
the M33 rotation curve. Columns are as follows: 1. parameter 
id; 2. best fit; 3., 4. mean and median from the marginalized 
likelihood; 5., 6. 68 and 95% confidence ranges. Note that we 
do not report constraints on ao for the motivations discussed 
in text. Finally, (Y*,ro,fXo) are in units of (Mq/Lq , kpc, km/s) 
respectively, while (77, ^C) are dimcnsionlcss numbers. 



Markov chain algorithm running two chains with 100000 
points each and ch ecking convergence throu gh the Gelmann - 
Rubin criterium ()Gelman fc Rubinl 1 19921 ) . The marginal- 
ized constraints on the model parameters are then obtained 
by the histograms of the merged chains after cutting out the 
burn in phase and thinning to avoid spurious correlations. 
The best fit parameters turns out to be : 



T. 



ao 



= 1.21 
0.0 , 



5.25 



K, = 26.27 



ro 



13.08 kpc 



(TO 



76.4 km/s 



giving an excellent fit to the data as shown in Fig.[7]and also 
confirmed by the value of the reduced \ 2 (x 2 /d.o.f. = 0.59). 
Although such a low \ 2 could also signal that the errors are 
overestimated (as it is possible given how they have been 
estimated) , the very good agreement between the theoretical 
and observed curve and the lack of any systematic trend of 
the residuals allow us confidently arguing that the model is 
indeed a good representation of the data. 

Table Q] provides the constraints on the model parame- 
ters from the rotation curve fit. Apparently, the confidence 
ranges are quite large for all the cases except the stellar M/L 
ratio. Indeed, this is a consequence of the severe degeneracy 
among the halo parameters which could be solved by adding 
some physically motivated priors. Table [1] does not report 
constraints on ao . Actually, our code gets stuck near ao = 
and starts wandering around this point because of a local 
minimum. We have tried to avoid this technical difficulty 
by running shorter chains for fixed values of ao finding out 
that the constraints on the other parameters are nearly un- 
affected and therefore concluding that it is not possible to 
constrain this quantity with the data at hand. Actually, such 
a result could be anticipated noting that all the quantities 
we have discussed up to now are essentially independent on 
this parameter unless one is able to probe the very inner 
regions (x « 0.01) which is not possible observationally. 
Since we are here interested in probing the viability of our 
model rather than constraining its parameters, we will not 
try to narrow down the above confidence ranges, but we plan 
to do this for a larger galaxy sample in a forthcoming work. 
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5 CONCLUSIONS 

Modelling galaxies is one of the classical topics of galactic 
dynamics. While empirical, observationally motivated and 
N-body simulation inspired, models are widely used, they 
nevertheless lack a derivation from an underlying theoretical 
approach and do not allow to understand the effect of the 
different phenomena shaping the mass density profile. As a 
classical alternative, one can rely on the analogy between 
stars in a galaxy and molecules in a gas to extend the meth- 
ods of statistical mechanics to the study of galactic systems. 
In this framework, we made a step forward investigating 
which density profile turns out as equilibrium configuration 
in the context of Tsallis nonextensive statistics. 

Under the assumption of spherical symmetry we demon- 
strated that power - law DFs can solve the generalized Boltz- 
mann equation and, when coupled to the Poisson equation, 
generate self-consistent models, which can be used to de- 
scribe dark matter haloes. Subsequently, we investigated in 
detail the properties of the underlying mass model showing 
that its attractive features provide an interesting alternative 
to empirical density laws. As a test case, we have also fitted 
with remarkable success the M33 rotation curve, providing a 
significant observational support to the theoretical context. 

A key ingredient within our derivation was the assump- 
tion of a parameterized expression for the velocity disper- 
sion term entering the power -law DF. Although the result- 
ing model is well behaved both theoretically and observa- 
tionally, it is nevertheless worth noticing that using Eq. (|14[) 
rather than another functional expression is actually only 
a matter of taste. Indeed, as far as we know, there are no 
theoretical motivations which can drive the choice of this 
function, wherefore it is only possible to introduce a heuris- 
tic approach and checking a posteriori the viability of a given 
cr(r) profile. Alternatively, one could assume an expression 
for q(r) and then solve Eqs, (|12|l - (|13|l to get both p(r) and 
a(r) but the choice of a functional form for q(r) is difficult 
since we can only argue that q(r) should approach the value 
of 1 asymptotically. The question of the central value of q(r) 
and how the asymptotic regime is approached is far to be 
understood, implying that this path seems to be affected by 
a large degree of arbitrariness. 

Regarding the uncertainty which approach to follow, the 
present proposal must be viewed as a first reasonable step re- 
quiring furhter investigations of the observational viability. 
Indeed, the successful fit of the M33 rotation curve appears 
as strong evidence in favour of the nonextensive approach, 
but is just a first step. In continuation it has to be checked 
whether this result is confirmed by fitting a larger sample 
of galaxies, spanning a large mass and luminosity range. On 
the other hand, rotation curves are hardly measured in el- 
liptical galaxies implying that for these systems one has to 
fit the velocity dispersion profile. Using both stars and, e.g., 
planetary nebulae allows to trace the line of sight velocity 
dispersion up to large radii, thus providing the possibility 
of both, checking whether the nonextensive model applies 
also for theses systems and constraining the parameters in- 
volved. Finally, proceeding to higher redshift, the combined 
fit to projected mass and aperture velocity dispersion in lens 
systems (such as the SLACS sample) should test not only 
the model, but offer also hints on the evolution of the un- 
derlying parameters. 



If all the above tests can be performed successfully, we 
would have not only a dark halo model inspired by funda- 
mental principles, but also a new recipe to construct, in a 
self consistent way, dynamical models for galaxies in view of 
further steps towards an understanding of the nature of one 
of the main components in our cosmos. 
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